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ABSTRACT 

Granular materials of different sizes are present on the surface of several atmosphere- 
less Solar System bodies. The phenomena related to granular materials have been 
studied in the framework of the discipline called Granular Physics; that has been 
studied experimentally in the laboratory and, in the last decades, by performing nu- 
merical simulations. The Discrete Element Method simulates the mechanical behavior 
of a media formed by a set of particles which interact through their contact points. 

The difficulty in reproducing vacuum and low-gravity environments makes numer- 
ical simulations the most promising technique in the study of granular media under 
these conditions. 

In this work, relevant processes in minor bodies of the Solar System are studied 
using the Discrete Element Method. Results of simulations of size segregation in low- 
gravity environments in the cases of the asteroids Eros and Itokawa are presented. The 
segregation of particles with different densities was analysed, in particular, the case of 
comet P/Hartley 2. The surface shaking in these different gravity environments could 
produce the ejection of particles from the surface at very low relative velocities. The 
shaking causing the above processes is due to: impacts, explosions like the release of 
energy by the liberation of internal stresses or the re accommodation of material. Sim- 
ulations of the passage of impact-induced seismic waves through a granular medium 
were also performed. 

We present several applications of the Discrete Element Methods for the study of 
the physical evolution of agglomerates of rocks under low-gravity environments. 

Key words: minor planets, asteroids: general - comets: general - methods: numerical 



1 INTRODUCTION 

Granular materials of different sizes are present on the sur- 
face of several atmosphere- less Solar System bodies. The 
presence of very fine particles on the surface of the Moon, the 
so-called regolith, was confirmed by the Apollo astronauts. 
From polarimetric observations and phase angle curves, it 
is possible to indirectly infer the presence of fine particles 
on the surface of asteroids and planetary satellites. More re- 
cently, the visit of spacecraft to several asteroids and comets 
has provided us with close pictures of the surface, where par- 
ticles of a wide size range from cm to hundreds of meters 
have been directly observed. The presence of even finer par- 
ticles on the visited bodies can also be inferred from image 
analysis. 

It has been proposed that several typical processes of 
granular materials, such as the size segregation of boulders 
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on Itokawa, th e displacement of boulders on Eros, among 
others (see e.g. lAsphaud l|2007h and references therein), can 
explain some features observed on the surfaces of these bod- 
ies. The conditions at the surface and the interior of these 
small Solar System bodies are very different compared to 
the conditions on the Earth's surface. Below we point out 
some of these differences: 



• while on the Earth's surface the acceleration of gravity 
is 9.8 m/s^ with minor variations, on the surface of elon- 
gated fcm-size asteroid is on the order of 10"^ to 10"* m/s^, 
with typical variations of a factor of 2 

• the presence of an atmosphere or any other fiuid media 
plays an important role on the evo lution of grains, partic- 
ularly in the small ones l|Pak et al.l (|l995i )). Under vacuum 
conditions in space, this effect does not occur. 

• In vacuum and low gravity conditions, other forces 
might play a role comparable to that of gravity, e.g. van der 
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Waals forces (jScheered l|2010l )). although these forces are not 
considered in our present approach. 

The phenomena related to granular material have been 
studied in the discipline called Granular Physics. Granular 
media are formed by a set of macroscopic objects (grains) 
which interact through temporal or permanent contacts. 
The range of materials studied by Granular Physics is very 
broad; rocks, sands, talc, natural and artificial powders, 
pills, etc. 

Granular materials show a variety of behaviours under 
different circumstances: when excited (fluidised), they often 
resemble a liquid, as is the case of grains flowing through 
pipes; or they may behave like a solid, like in a dune or a 
heap of sand. 

These processes have been studied experimentally in 
the laboratory, and, in the last decades, by numerical anal- 
ysis. The numerical simulation of the evolution of granular 
materials has been done recently with the Discrete Element 
Method (DEM) . DEM is a family of numerical methods for 
computing the motion of a large number of particles such as 
molecules or grains under given physical laws. DEMs simu- 
late the mechanical behavior in a media formed by a set of 
particles which interact through their contact points. 

Low-gravity environments in space are difficult to re- 
produce in a ground-based laboratory; especially if one is 
interested in keeping a stable value of the acceleration of 
gravity on the order of 10~^ to 10~* m/s^ for several hours, 
since under these low-gravity conditions the dynamical pro- 
cesses are much slower than on Earth. Parabolic flights are 
not suitable for these experiments, since it is not possible 
to attain a stable value during the free-fall flight. For labo- 
ratory experiments, we are then left with experiences to be 
carried on board space stations. 

Therefore, numerical simulation is the most promising 
technique to study the phenomena affecting granular mate- 
rial in vacuum and low-gravity environments. 

The rest of the article is organised as follows, fn Section 
[5] we describe the implementation of the Discrete Element 
Methods used in our simulations. In Section |3] we present 
the results of simulations of the process of size segregation 
in low-gravity environments, the so-called Brazil nut effect, 
in the cases of Eros, Itokawa and P/Hartley 2. In Section |4l 
the segregation of particles with different densities is anal- 
ysed, with the application to the case of P/Hartley 2. The 
surface shaking in these different gravity environments could 
produce the ejection of particles from the surface at very low 
relative velocities; this issue is discussed in Section [5] The 
shaking that causes the above processes is due to impacts 
or explosions like the release of energy by the liberation of 
internal stresses or the re accommodation of material. Al- 
though DEM methods are not suitable to reproduce the im- 
pact event, we are able to make simulations of the passage 
of impact-induced seismic waves through a granular media; 
these experiments are shown in Section [6] The conclusions 
and the applications of these results are discussed in Section 

m 



2 DISCRETE ELEMENT METHODS 

DEM are a set of numerical calculations based on statisti- 
cal mechanic methods. This technique is used to describe 



the movements of a large amount of particles which are sub- 
jected to certain physical interactions. 

DEMs present the following basic properties that gen- 
erally define this class of numerical algorithms: 

• The quantities are calculated at points fixed to the ma- 
terial. DEM is a case of a Lagrangian numerical method. 

• The particles can move independently or they can have 
bounds, and they interact in the contact zones through dif- 
ferent types of physical laws. 

• Each particle is considered a rigid body, subject to the 
laws of rigid body mechanics. 

The forces acting on a particle are calculated from the in- 
teraction of this particle with its nearest neighbors, i.e. the 
particles it touches. Several types of forces are usually con- 
sidered in the literature; e.g free elastic forces, bonded elastic 
forces, frictional forces, viscoelastic forces, interaction of the 
particles with other objects, such as walls and mesh objects 
acting as boundary conditions, global force fields (i.e. grav- 
ity), velocity dependent damping, etc. 

The main drawback of the method is the computational 
cost of computing the interacting forces for each particle at 
each time step. A simple all-to-all approach would require to 
perform 0{N{N — l)/2) operations per time step, where A'^ 
is the number of particles in the simulation. Several efficient 
methods to reduce the number of pairs to compute have 
been implemented; e.g. the Verlet lists method, the link cells 
algorithm, and the lattice algorithm. Another problem for 
the simulation is the length of the time step, which should 
be much less than the duration of the collisions, typically 
1/10 to 1/20 of collisions duration. Based on the Hertzian 
elastic contact theory, the duration of contact (r) can be 
expressed as: 



r = 5.84 



E 



(1) 



jWada et al.l (|2006l '). after iTimoshenko and Goodieij 
(ll97ol )'). where o is the grain density, nu is the Poisson ratio, 
E is the material strength, r is the radius of the particle, and 
V is the coUisional velocity. 

In Figure [1] we plot the previous estimate of the dura- 
tion of collision as a function of the coUisional velocity for 
particles with r = 0.1, 1 and 10 m. The other parameters 
are assumed as follows: p — 2000 gr/crn^, v = 0.17, and 
E = 100 Gpa. We are interested in the processes that oc- 
cur on the surface of the small Solar System bodies, where 
the interactions among the boulders occur at velocities com- 
parable to the escape velocity on their surface. The upper 
i-axis indicates the radius of the body, while the lower one 
shows the corresponding escape velocity (assuming a con- 
stant density of p = 3000 gr/cm^). For fcm-size asteroids 
and m-size boulders, the collisions typically last for a few 
10~''s, therefore, the time step required to correctly simu- 
late the collision would be ~ 10~'*s. 



2.1 Viscolelastic spheres with friction 

The contact force between two spherical particles can be de- 
composed in two vectors (Figure [5}: the normal force, along 
the direction that joins the centres of the interacting par- 
ticles; and the tangential force, perpendicular to this line. 
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Figure 1. Estimate of the duration of collision as a function of 
the collisional velocity for particles of r = 0.1, 1 and 10 m 
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Figure 2. Scheme of the contact forces between two spherical 
particles 



Naming i and j the two interacting particles, the total force 
Fij can then be expressed as: 



F, 



if i/'jj > 
otherwise 



(2) 



where F^j is the normal force and F^j the tangential 
one. ^pij is the deformation given by: 



Ri + Rj 



(3) 



where Ri and Rj are the radius of the particle i and j, 
respectively, and rt and rj' are the position vectors. 

Several models have been used for the normal and tan- 
gential forces in the literature. Among the most used ones 
is the damped dash-pot, also known as the Kelvin- Voigt 
model. Instead of using this model in our simulation, we use 
an ex tension of an elastic-spheres one developed by iHert j 
()l882l ). since it is a more realistic representation of two col- 
liding particles. 

The normal interaction force between two elastic 
spheres -Fj"'"^' was inferred by Hertz as a function of the 
deformation ip: 



3(1 



(4) 



where Y is the Young modulus and v is the Poisson 
ratio. The effective radius Ref / is given by the expression: 



R, 



eff 



1^1 
Ri Ri 



(5) 



A viscoelastic interaction between the particles can be 
modelled by including a dissipation factor in eq. ID The vis- 
coelastic normal forces F"'^'^ then become: 



jpn;ve 



2Yy/R, 



3(1 



(6) 



where ^4 is a dissipative constant and dil^/dt is the time 
derivative of the deformation. 

Considering the previous expression for the normal force 
could lead to unrealistic results, since it does not take into 
account the fact that the particle s do not overlap, but they 
become deformed (Poschel and Schwageij (2005)). During 
the compression phase and most of the decompression phase, 
the term {^j^^^ + Ayftp^'^ in eq. [6] is positive, leading to a 
repulsive (positive) normal force. However, at a certain stage 
of the decompression, the deformation ip could still be pos- 
itive, but the second term could be negative, which would 
lead to a negative (attractive) force. This is an unrealistic 
situation, since there are no attractive forces during the col- 
lision of two particles. The problem arises when the centres 
of the particles separate too fast from one another to allow 
their surfaces to keep in touch while recovering their shape. 
In order to overcome this problem, for the condition tj) > Q 
in eq.[5J we use the following expression for the normal force: 



0, 



3(1 -;/2) 



Following the model by ICundall and StrakI l| 19791 ') for 
the tangential force (-F*), when two particles first touch, 
a shear spring is created at the contact point. The static 
friction is then modelled as a spring acting in a direction 
tangential to the contact plane. The particles start sliding 
with the shear spring resisting the motion. When the shear 
force exceeds the normal force multiplied by the friction co- 
efficient, dynamic sliding starts. We limit the shear force by 
Coulomb's friction law; i.e. |F* < /i_F"''"'^|. The expression 
for the tangential force then becomes: 



F* = -sign{vl.^i) min{||K?|| , mII-P""'""!!} 



(8) 



where the first term inside the curly brackets corre- 
sponds to the static friction, and the second one is the dy- 
namic friction, k is a constant, ? is the elongation of the 
spring, and jj, is the dynamic friction parameter. 



2.2 ESyS-particle 

For the DEM simulations we d eveloped a versio n 
of the ESyS-particle package (|Abe et al\ (|2004 ): 
https://launchpad.net/esys-particle) adapted to our needs. 
ESyS-particle is an Open Source software for particle-based 
numerical modeling, designed for execution on parallel 
supercomputers, clusters or multi-core computers running a 
Linux-based operating system. The C-|— |- simulation engine 
implements a spatial domain decomposition for parallel 
programming via the Message Passing Interface (MPI). A 
Python wrapper API provides flexibility in the design of 
numerical models, specification of modeling parameters and 
contact logic, and analysis of simulation data. 

The separation of the pre-processing, simulation and 
post-processing tools facilitates the ESyS-particle develop- 
ment and maintenance. The setup of the model geometry is 
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given by scripts, since the whole package is script driven (no 
interactive GUI is provided by ESyS). 

The particles can be either rotational or non-rotational 
spheres. The material properties of the simulated solids can 
be elastic, viscoelastic, brittle or frictional. Particles can 
be bonded to other particles in order to simulate break- 
able material. It is possible to implement triangular meshes 
for specifying boundary conditions and walls. The package 
also includes a variety of particle-particle and particle-wall 
interaction laws; such as linear elastic repulsion between un- 
bounded contacting particles, linear elastic bonds between 
bonded particle pairs, non-rotational and rotational fric- 
tional interactions between unbounded particles, rotational 
bonds implementing torsion and bending stiffness and nor- 
mal and shear stiffness. Boundary conditions and walls can 
move according to pre-defined laws. 

The DEM implementation in ESyS-particle employs the 
explicit integration approach, i.e. the calculation of the state 
of the model at a given time only considers data from the 
state of the model at earlier times. Although the explicit 
approach requires shorter time steps, it is easier to develop a 
parallel version for execution in high performance computing 
inf rast ruct ures . 

ESyS-particle has shown good scaling performance 
when using additional computing elements (processor cores) , 
if at least ~ 5000 particles are processed by each core. Oth- 
erwise, when a lower number of particles is handled by each 
core, the impact of the overhead by the communications be- 
tween processes reduces the computational efficiency of the 
application. As long as the problem size is scaled with the 
number of cores, the scalability is close to linear. Therefore, 
large amounts of particles, typically a few million, are pos- 
sible to model. 

For the analysis of the results, ESyS-particle can for- 
mat the output to be used in 3D visualisation platforms like 
VTK and POV-Ray. In particular, we use the software Par- 
aview, based on VTK and developed by Kitware Inc. and 
Sandia National Labs (EEUU), which offers good quality in 
3D graphics and allows us to implement several visualisation 
filters to the data. 

ESyS-particle has been employed to simulate earth- 
quake nucleation, comminution in shear cells, silo flow, rock 
fragmentation, and fault gouge evolution, to name but a few 
applications. Just to give a few references, we mention exam- 
ples in fr acture mechani c s (Schopfer et al.. (.200&)). fau lt me- 
chanics l|Abe and Maiij (|2005iriMair and Abel (|2008l ) ) . and 
fault rupture propagation ( Abe et al. ( 20061 )'). 

For our simulations, we have implemented the Hertzian 
viscoelastic interaction model with and without friction into 
the ESyS-particle package, according to the eqs. [7] and 
[8] Several modifications were necessary to implement ei- 
ther in the C-|— I- code as well as in the Python interface 
(|Heredia and Richeril (|2009h '). 



2.3 Tests 

In order to test the code and to set the values of the rele- 
vant physical and simulation parameters, we choose a few 
problems for which there exists an analytical solution or we 
can compare the output with experiments. 



2.3.1 Test case 1: a direct collision of two equal spheres 

We consider the case of two equal viscoelastic spheres: one 
starts at rest and the other one approaches from the nega- 
tive a;-direction at a given speed along the line joining the 
particle centres. Friction is not considered, since the colli- 
sion between the particles is normal. The aim of this test is 
to study the viscoelastic collision. 

The coefficient of restitution can be used to characterise 
the change of relative velocity of inelastically colliding parti- 
cles. Let us note and the velocities before the collision 
of particles 1 and 2, respectively; and and V2 the veloc- 
ities immediately after. When the relative velocity is along 
the line joining the particle centres, we note v = \vt — vt\ 
and v' — [iPi The coefficient of restitution e is then 

calculated as: 



(9) 



In general, this coefficient depends not only on the im- 
pact velocity, but also on material properties. 

Because of their deformation, particles lose contact 
slightly before the distance of the centres between the 
sphere s reaches the sum of the radii. ISchwager and Poschell 
(2008") present an analytical estimate of the coefficient of 
restitution which takes into account this fact. The computa- 
tion of e is then presented as a divergent series of the dimen- 



sionless parameter Pv^^^\ where P 



3 pA 



2 m, 



5 = 



-I- 



— . The material parameters Y, v and A 



are defined above. R-i, mi, m2 are the radius and mass 
of particle 1 and 2, respectively. 

We run simulations of two colliding particles with the 
following combination of parameters: Y — {10^,10^°} Pa, 
A = {10~*, 10""^} V — 0.3; for a set of initial relative 
velocities v — {0.1,0.5,1,5,10} m/s. The particles have a 
radii of Im and a density of 3000 kg/m? . The time step of 
the integration is 10~^ s. 

The coefficient of restitution for the numerical simula- 
tions is presented in Figure |3^ as a function of Pv^^^. The 
symbols correspond to different combinations of parameters: 
circle - Y = 10^ A = 10"*; down triangle - Y = 10^ 



A — 10 ^; square 



Y = 10' 



10 ; up triangle 



Y = 10"', A = 10-^ {Y in Pa and A in s"'). The analyt- 
ic al estimates are compu t ed wi th Maple's codes presented 
in I Schwager and Poschell (|2008l ) , where the expansions are 
up to 40th order. The ratio between the numerical and an- 
alytical estimate is presented in Figure [Sjs. We find a good 
agreement between the two estimates up to values of Pv^^^ 
closer to 1. The discrepancy is due to the cut-off of the higher 
terms. 

Several laboratory experiments have been conducted to 
esti mate the coefficie n t of restitution of ro ck materials (see 
e.g. llmre efai] (|2008h . lDurda et al\ (|201ll ')'). For impact ve- 
locities in the range 1 — 2 m/s, values of e ~ 0.8 — 0.9 
have been obtained. Looking back at Figure [3^, we observe 
that this range of values of e are obtained for the following 
set of material parameters: Y = lO'^" Pa, A — s~^, 
= 0.3. Therefore, we will choose these parameter values 
for our numerical simulations of colliding rocky spheres. 
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Figure 3. a) The coefficient of restitution for the numerical sim- 
ulations of the collision of two viscoelastic spheres as a function of 
fiv^/^. The symbols correspond to different combinations of pa- 
rameters: circle — Y = 10^, A = 10~*; down triangle — Y = 10^, 
A = 10-3; square -Y = 10^°, A = 10"*; up triangle -Y = 10^", 
A = 10 -3 (y in Pa and A in s'^). b) The ratio between the 
numerical and analytical estimate. 



2.3.2 Test case 2: a grazing collision between two spheres 

In this case we consider two equal viscoelastic spheres: one 
starts at rest and the other one approaches from the negative 
x-direction at a given speed; but, in contrast to the previ- 
ous case, the distance between the y- values of the particles 
centres is slightly less than the sum of the radius. We then 
have a grazing collision. The aim of this test is to compare 
the results of the viscoelastic interaction with and without 
friction. We run simulations of two colliding particles with 
the following set of parameters: Y — 10^" Pa, A = IQ-^ s~^, 
V — 0.3, i?i = -R2 = 1 Ti, and a density of p = 3000 kgjrrC' . 
The friction parameters of eq. |H] are chosen as; k = 0.4, 
/i = 0.6. The distance between the centres in the j/-direction 
is (0.999-Ri -f -R2). We run simulations where particle 2 has 
initial velocities oiv = {IQ-'^, 0.01, 0.1, 1, 10} m/s. 

In Figure 3] we plot the ratio between the modulus of 
the particles relative velocity after exiting and before the 
interaction, as a function of the initial velocity. The star 
symbols correspond to the simulations without the friction 
interaction and the cross symbols to the ones with it. Due 
to the fact that the collision is almost grazing, the ratios 
are almost 1 for the simulations without friction, regardless 
of the initial velocity. For simulations with friction, as it is 
expected, the ratio decreases as the initial velocity decreases, 
because the friction interaction becomes more relevant for 
low velocities. 



2.3.3 Test case 3: a bouncing ball 

In ESyS-particle the interaction between a particle and a 
mesh wall can be linear, elastic or a linear elastic bond. 
Viscoelastic and frictional interactions of particles and walls 
are not yet implemented. Therefore, in order to simulate 
a frictional viscoelastic collision of a ball against a fixed 
wall, we have to glue balls to the wall with a linear elastic 
bond. The following test case consists on a free-falling ball 
impacting on an equal size ball that is bonded to the floor. 
The objective of this experiment is to test different time 
steps for the simulations. 

Y = 10^° Pa, 
m, and a density of 
The bonded particle has an elastic bond 



We use the following set of parameters 
A = 10-3 s-^ u = 0.3, Ri = R2 = 1 
p = 3000 kg/m^ 
with a modulus K = lO'' Pa. Particle 2 falls from a height of 




10 10 

Initial Velocity (m/s) 

Figure 4. The ratio between the modulus of the relative velocity 
between the particles after exiting and before the interaction as a 
function of the initial velocity. The symbols correspond to differ- 
ent contact force models: square - Hertzian viscoelastic spheres; 
circle - Hertzian viscoelastic spheres with friction. 
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Figure 5. a) The distance of the falling particle respect to the 
edge of the resting one (centre height minus Si?) as a function 
of time, b) The ratio of the previous values to the one for the 
smallest time step at each time. 



2.75m. In the first set of simulations we assume the Earth's 
surface gravity {g = 9.81 m/s"^). For this set, we use the 
following time steps in the simulations: = {6 x 10-**, 5 x 
10-*, 10-*, 10-^ 10-^ 5 X lO-''} s. 

The duration of the collisions is computed from the sim- 
ulations as the interval of time while the deformation param- 
eter defined in eq. [3]is greater than 0. As mentioned above 
this interval is slightly larger than the time the balls are in 
contact, but it is good enough for the purpose of having an 
order of magnitude estimate of it. For the previous set of 
parameters, the duration of the collision is ~ 0.003 s. 

In Figure [5^, we plot the distance of the falling par- 
ticle respect to the edge of the resting one (centre height 
minus 3-R) as a function of time. In Figure [SJs, we plot the 
ratio of the previous values to the one for the smallest time 
step at each time. We find that for time steps dt ^ 10-^ s, 
there is a very good agreement between the simulations. For 
longer time steps the bouncing ball presents an implausible 
behavior. 

The coefficient of restitution can be computed as the 
ratio between the velocity at the iteration step just after 
the collision and at the step just before the collision (just 
after and before the deformation defined in eq. |3]is 1/; < 0). 
For time steps dt ^ 10-^ s there is a good agreement among 
the different estimates. We obtained a value of 0.593. 
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Therefore, for the previous set of parameters, we wiU 
use a time step of dt — 10~^ s for the simulations with 
Earth's gravity, since the coUision is covered with ~30 time 
steps and it is a good compromise between quahty of the 
resuhs and a longer time step. 

In another set of simulations we use very low surface 
gravity, similar to the one found on the surface of aster- 
oid Itokawa and comet P/Hartley 2; i.e. a rocky object of 
~ 500 m in diameter or an icy object of ~ 1 km in di- 
ameter. For this set, we use the following time steps in the 
simulations: = {5 x 10"*, 10"*, 10"^ lO"**} s. For the 
previous set of parameters, the duration of the collision is 
~ 0.01 s. For time steps dt 10~* s there is a good agree- 
ment among the different runs. The coefficient of restitution 
in these simulations is 0.721. For the simulations in this low- 
gravity environments we will use a time step dt — 10~* s, 
which corresponds to a collision lasting ~ 100 time steps. 

2. 3. 4 Test case 4: Newton's cradle 

A Newton's cradle is a device used to demonstrate the con- 
servation of linear momentum and energy via a series of 
swinging hard spheres. When one ball at the end is lifted 
and released, it knocks a second ball and this one the next 
until the last ball in the line is pushed upward. A typical 
Newton's cradle consists of a series of identically sized metal 
balls hanging by equal length strings from a metal frame so 
that they are just touching each other at rest. 

We simulate the Newton's cradle with four spheres 
aligned in the x-axis. We number the particles from right 
to left: #1 being the particle at the right extreme and #4 
the one at the left extreme. The x-axis increases to the right. 
Particle #1 has a negative initial velocity Vx = —10 m/s. 
Two types of simulation are run: Hertzian elastic and vis- 
coelastic spheres. We use the following set of parameters: 
Y = 10^° Pa, A = 10"^ s~\ = 0.3 (for the viscoelastic 
simulation). The radius of the spheres are R — 1 m, and a 
density oi p — 3000 kg/m? . 

In Figure [5] we present the time evolution of the fol- 
lowing parameters for each simulation: i) a;-position of each 
particle (#1 to #4); ii) x-velocity for each particle; iii) 
relative change of x-total momentum: {Momentum{t) — 
Momentum{t = 0)) / M omentum{t = 0); iv) relative change 
of total kinetic energy: {K.E{t)- K.E .{t = Q))/K.E.{t = 0). 
Figure [S] a) corresponds to the Herztian elastic (HE) simu- 
lation, and Figure |5] b) to the Herztian viscoelastic (HVE) 
one. 

Note that for the HE simulation particle #4 acquires 
almost the velocity of the initial impacting particle and lit- 
tle rebound is observed in the particles ^1 to #3. The linear 
momentum is conserved after the collision up to a relative 
precision < 10~^^, and the kinetic energy after the rebound 
is conserved up to a relative precision of 10~^^. In the HVE 
simulation, the particle #4 acquires 70% of the velocity of 
the initial impacting particle, and particle #3 acquires 25%. 
No rebound is observed and all the particles move to the 
left. The final velocities increase from right to left. The lin- 
ear momentum is also conserved after the collision up to a 
relative precision < 10~^^ (down to the last output digit). 
The kinetic energy after the rebound is not conserved ~ 50% 
of the initial kinetic energy is spent on the damping of the 
viscosity interaction. 
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Figure 6. a) Results of the Herztian elastic (HE) simulations of 
the Newton's cradle: i) x-position of each particle (#1 to #4); ii) 
x-velocity for each particle; iii) relative change of x-total momen- 
tum; iv) relative change of the total kinetic energy, b) Similar set 
of plots for the Herztian viscoelastic (HVE) simulations of the 
Newton's cradle. 



3 SIZE SEGREGATION IN LOW-GRAVITY 
ENVIRONMENTS: THE BRAZIL NUT 
EFFECT 

3.1 The shaking or knocking procedure 

Consider a recipient with one large ball on the bottom and a 
number of smaller ones on top of it. All the balls have similar 
densities. After shaking the recipient for a while, the large 
ball rises to the top and the small ones s i nk to the bottom 
(|Rosato et all (|l987l ). lKnight et al.l (fl99l ). lKudro"iiil l|2004h ). 
This is the so called Brazil nut effect (BNE), because it can 
be easily seen when one mixes nuts of different sizes in a can; 
the large Brazil nuts rise to the top of the can. Unless there 
is a large difference in the density of the balls, a mixture of 
different particles will segregate by size when shaken. 
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Figure 7. The floor is vcrticaUy displaced at a certain speed 
(ufloor) for short interval (rftshafce)i according to a staircase- 
like function. 

The BN E has b een attributed to the foUowing processes 
ijHong et al\ (j200ll )): i) the percolation effect, where the 
small er ones pass through the holes created by the larger 
ones l|jullien and MeakinI ([l99^)); ii) geometrical reorgani- 
sation, through which small parti cles readily fil l small open- 
ings below the large particles ijRosato et alJ ()l987l )): iii) 
global convection which brings the large par ticles up but 
does n ot allow for reentry in the downstream dKnight et alJ 
l| 19931 ) 1: iv) due to its larger kinetic energy, the large par- 
ticle still foll ows a baUistic upraise, penetra ting by inertia 
into the bed (|Nahmad-Mohnari etall HqoI)). 

While size ratio is a dominant factor, particle-specific 
properties such as density, inelasticity and friction can also 
play important roles . 

IWilliam^ (|l963l ') performed a model experiment with a 
single large particle (intruder) and a set of smaller beads 
inside a rectangular container. When the container was vi- 
brated appropriately, the intruder would always rise and 
reach a height in the bed that depends on vibration strength. 

In order to simulate this effect under different gravity 
conditions, we run simulations of a 3D box with many small 
particles and one big p article at t he bo t tom, the s o-call ed 
intruder model system (jWiUiamd (|l963h . lKudrdiil (|2004 )). 
On the fioor we glue one row of small particles with a linear 
elastic bond. The box is subjected to a given surface gravity. 

We run simulations under several gravity conditions: 
the surface of the Earth, Moon, Ceres, Eros and a very-low 
gravity environment like the surface of asteroid Itokawa or 
comet P/Hartley 2. The parameters for the simulations are 
summarised in Table [T] The physical and elastic parameters 
of the particles are similar to the ones used in the previous 
tests: Y = 10^° Pa, A = 10"^ u = 0.3, k = 0.4, fi = 0.6, 
K = 10^ Pa, p = 3000 kg/m^. 

The fioor is vertically displaced at a certain speed 
(vfioor) for a short interval (dtshake), according to a 
staircase-like function like the one presented in Fig [T] The 
process is repeated every given number of seconds (Atrep), 
depending on the settling time given by the surface gravity. 
We have chosen this vibration scheme instead of the fre- 
quently used sinusoidal oscillation of the fioor, because we 
are interested in the effects of a sudden shock coming from 
below. This shock could arise from the translation of the im- 
pulse generated by an impact in a far region. We refer this 
vibration scheme as a shaking or knocking procedure. 

In order to prepare the initial conditions for the simu- 



lations of the BNE, we run a set of simulations where the 
particles start at a certain height over the surface and they 
free fall. The fioor is slightly shaken at the beginning of 
these preliminary simulations in order to obtain a random 
settling of the particles. After finishing the shaking and let- 
ting the particles settle down, we use the positions at the 
end of the runs as the initial conditions for the set of BNE 
simulations. We must run different preliminary simulations 
for each gravity environment. 

In the BNE simulations, the floor's velocity is linearly 
increased from up to the final value v floor, which is reached 
after 20 jumps. We note that the shaking procedure is pa- 
rameterized with the floor's velocity. 

The 3D box is constructed with elastic mesh walls. The 
box has a base of 6 x 6 m and a height of 150 m. A set of 
12 X 12 m small balls of radius Ri = 0.25 m are glued to the 
floor. The big ball has a radius R2 = 0.75 m, and on top of it, 
there are 1000 small balls with a normal distribution of radii 
(mean radius Ri — 0.25 m, standard deviation a = 0.01 m). 
We use the same box for all the simulations. 

The size range of the balls are selected in correspon- 
dence with the boulders size observed on the surface of as- 
teroid Itokawa and Eros. 



3.2 Earth 

We run simulations with the following set of floor veloci- 
ties: V floor = {0.3,1,3,5,10} m/s. Snapshots at start and 
after 100 shakes (100 sec. of simulated time) are presented 
in Figure ID The snapshots correspond to the simulation 
with floor's velocity v floor ~ 5 m/s. In the supplementary 
material we include movies with the complete simulation 
(moviel with all the spheres drawn and movie2 with the 
small spheres erased). 

In Figure [5] we present the evolution of the big ball's 
height as a function of the number of shakes for the differ- 
ent floor velocities. The thick black line marks the height 
of a box enclosing the 1000 small particles with a random 
close packing. Ra ndom close packin g has a max;imum poros- 
ity of P = 0.64 (| Jaeger and Nagell ^199^ )). The volume of 
the enclosing box is calculated as the sum of the volume 
of the 1000 small particles divided by the porosity; i.e.: 
V = 1000{^iRf)/P = 102 For a box with a 6 x 6 m 
base, we obtain a height of the enclosing box of 2.84 m. The 
thick black line is drawn at this height. 

For the two lowest velocities {v floor ~ {0.3, 1} m/s) 
the big ball stays at the bottom, for the two largest ones 
{v floor = {5, 10} m/s) it rises to the top, and for the inter- 
mediate one {v floor ~ S m/s) it starts rising but does not 
reach the top at the end of the simulation. 

When the floor's displacement velocity is below ~ 
3 m/s, the Brazil nut effect does not occur. Above this 
threshold, the time required by the big ball to reach the top 
decreases for increasing floor velocities. Note that there is 
a sharp decrease in the rising time for small changes in the 
floor's velocity (from 3 to 5 m/s). For large displacement 
velocities, the balls on the top, including the big one that 
is 27 times more massive than the small ones, can be lifted 
at considerable heights, as it is seen in the large excursions 
made by the big ball for v floor ~ 10 m/s. 
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Table 1. Parameters for the simulations of the BNE under different gravity environments. 



Parameter 




Earth 


Moon 


Ceres 


Eros 


Low-gravity 
Itokawa & 
P /Hartley 2 


Surface gravity g (m/s^) 




9.81 


1.62 


0.27 


5.9 X 10-3 


10-4 


Escape velocity Vesc (m/s) 




11.2 X 10^ 


2.38 X 10^ 


510 


10 


0.17 


Floor's velocity Vfi^or (jn/s) 




0.3 - 10 


0.1 - 3 


0.03 - 1 


0.01 - 0.3 


0.003 - 0.1 


Duration of displacement dtgi^ake 




0.1 


0.1 


0.1 


0.1 


0.1 


Time between displacements Ar-ep 


(s) 


2 


5 


5 


15 


15 



lii 




a) Moon 



b) Ceres 



Figure 8. Snapshots at start and after 50 shakes (100 sec. of sim- 
ulated time) for the simulation under Earth's gravity. The snap- 
shots correspond to the simulation with floor's velocity Vfi^or = 
5 m/s. See the movies in the supplementary material. 




20 30 

Number of shakes 

Figure 9. The evolution of the big ball's height as a function 
of the number of shakes for different floor velocities (f j-joor = 
{0.3, 1, 3, 5, 10} m/s) under Earth's gravity. Note that the floor's 
velocity is used as the varying parameter in the shaking process. 
A black line is drawn at a height of 2.84 m, which is the height 
of a compact enclosing box (see text). 




Number ot shakes 



150 200 



c) Eros 
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Figure 10. The evolution of the big ball's height as a function of 
the number of shakes for different floor velocities under different 
gravity environments: a) Moon; b) Ceres; c) Eros; d) Itokawa. 
The legends correspond to the floor velocities (f/ioor)- 



3.3 Comparison with other gravity environments 

Similar simulations were run for otlier gravity environments, 
like tlie surface of the Moon, Ceres, Eros and a very-low 
gravity environment like the surface of asteroid Itokawa 
or comet P/Hartley 2. The simulation parameters are pre- 
sented in Table [T] 

For the simulation under the very-low gravity environ- 
ment, we present a movie of 4500 sec. of simulated time 
(300 shakes) in the supplementary material. The movie cor- 
responds to the simulation with floor's velocity v floor = 
0.05 m/s (movieS with all the spheres drawn and movie4 
with the small spheres erased). 

Figure [To] presents the evolution of the big ball's height 
as a function of the number of shakes for the different floor 
velocities and the different gravity environments: a) Moon, 
b) Ceres, c) Eros, d) Itokawa. 

As in the cases of the simulations in Earth's gravity, in 
all the different gravity environments we can find a thresh- 
old for the floor's velocity, below which the Brazil nut effect 
does not occur. From the previous plots, we get a rough es- 
timate of these thresholds. In Figure [TT] we plot the velocity 
thresholds as a function of the surface gravity in a log-log 
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Figure 11. Comparison of the floor's velocity tlircsliold for the 
different gravity environments. The Brazil nut effect does not oc- 
cur if the floor's velocity is below the threshold. The floor's ve- 
locity thresholds are presented as small circles. Note that the 
thresholds are not precisely estimate, because they are computed 
from the plots in Figure[9land ll0l a-d. A straight line in the log-log 
space is fitted to the data points. The up triangles represent the 
escape velocity for the given surface gravity. The escape velocity 
for the largest objects are out of the plots. 



scale. A straight line in the log-log space is a good fit to the 
data points: 

logio vthr^ [ m/s] = 0.42 logio 9 [ "^/s^] + 0.05 (10) 

We conclude that the Brazil nut effect is effective in a 
wide range of gravity environments, expanding 5 orders of 
magnitude on surface gravity. 

In Figure 1111 we plot the escape velocity for the given 
surface gravity. Note that the floor's velocity thresholds ap- 
proach the escape velocity for the low gravity environments. 
For example, in the case of Itokawa, the escape velocity is 
fesc = 0.17 m/s, while the estimated floor's velocity thresh- 
old is Vfioor = 0.015 m/s. This point is revisited in Section 

El 



4 DENSITY SEGREGATION IN 

LOW-GRAVITY ENVIRONMENTS 

As mentioned above, other particle-specific properties can 
affect the segregation process. In particular, the effects of 
density have been studied the most. For ratios of the den- 
sity of the large to the small particles much larger than 1 
(denser large particles), the segregation effect could be re- 
versed, and the large particles would sinks to the bottom, 
producing the so-called Reve r se Brazil Nut Effe ct (RBNE) 
ijShinbrot and Muzzid l|l998t ). lHong et al\ (|200ll )). 

However, for particles of similar si z es bu t different 
densit ies, both labor atory (iMobius et al.l (|200ll ). IShi et all 
l|2007l )') or numerical l|Liml (|2010l )) experiments have shown 
that the lighter particles tend to rise and form a pure layer on 
the top of the system, while the heavier particles and some 
of the lighter ones stay at the bottom and form a mixed 
layer. In the Solar System, we might encounter bodies with 
such a mixture of heavy and light particles. Cometary nuclei 
are believed to be formed of a mix of icy and rocky material. 
However, the intimacy of this mixture is still unknown, with 





Figure 12. Snapshots of the initial and final state (after 1300 
shakes) for a simulation under the low-gravity environment and 
a fioor velocity oivfi^or = 0.05 m/s. (see movie5 and movie6 in 
the supplementary material. 



two possible scenarios: 1) every particle is made of a mix- 
ture of ice and dust, and 2) there exist some particles mainly 
formed by icy material and some others mainly formed by 
rocky constituents that are mixed together. 

We shall investigate the behavior of a mixture of light 
and heavy particles under different gravity environments. 

For the simulations we create a 3D box similar to the 
previous one, with a 6 x 6 m base and a height of 150 m. 
The box is constructed with elastic mesh walls. On the fioor 
we glue a set of 12 x 12 small balls of radius R\ = 0.25 m 
and density p = 2000 kg/m'^. There are 500 light balls with 
a normal distribution of radii (mean radius Ri — 0.25 m, 
standard deviation a — 0.01 m) and density p — 500 kg/m!^ . 
On top of them, there are 500 heavy balls with the same dis- 
tribution of radii and density p = 2000 kg/m^. At the be- 
ginning of the simulations the balls are placed sparsely, the 
light balls at the bottom and the heavy ones on top. They 
free fall and settle down before starting the fioor shaking. 

Elastic parameters of the particles are the same for both 
types of particles and similar to the ones used in the previous 
tests for aU the particles: Y = 10^° Pa, A = 10"^ s~\ 
u = 0.3, K = 0.4, p. = 0.6, K = 10^ Pa. 

The fioor is displaced with a staircase function in a sim- 
ilar way as in the previous set of simulations. 

Two gravity environments were tested: the Earth's sur- 
face gravity and a very-low gravity environment like the sur- 
face of comet P/Hartley 2. 

Figure[T2]presents snapshots of the initial and final state 
(after 1300 shakes) for a simulation under the low-gravity 
environment and a floor velocity of vfioor = 0.05 m/s. In 
the supplementary material we include movies with the com- 
plete simulations (movie5 corresponds to the simulation un- 
der Earth's gravity and v floor ~ 3 m/s; movie6 corresponds 
to the simulation under low gravity and v floor = 0.05 m/s. 
Note that in these movies the camera moves with the floor, 
therefore it seems that the floor is always located in the same 
position, but it really is moving with the staircase function 
described above). 

At every snapshot, we compute the median height of the 
light and heavy particles, respectively. These median heights 
are plotted as a function of the number of shakes in Figure 
I13l al for the Earth's gravity simulations, and b) for the low- 
gravity ones. For each simulation there are two lines: the one 
that starts on top corresponds to the heavy particles and the 
one that starts at the bottom to the light ones. In the Earth 
environment simulations, the lines do not cross for the two 
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Figure 13. The median height of the hght and heavy particles 
are plotted as a function of the number of shakes for different floor 
velocities and under different gravity environments: a) Earth, b) 
low-gravity like P /Hartley 2. The legends correspond to the floor 
velocities (t^/iooi )- For each simulation there are two lines: the 
one that starts on top corresponds to the medium height of the 
heavy particles (labeled with h) and the one that starts at the 
bottom to the medium height of the light ones (labeled with I) 



lowest floor velocities: u/joor = {li3} m/s; therefore, the 
particles do not overturn the initial segregation. Though, for 
Vfioor = 3 m/s, the lines start to approach. However, for 
the highest floor velocities, i.e. Vfioor = 5 m/s, the lines 
cross at an early stage of the simulation after which they 
remain almost parallel. Most of the hght particles move to 
the top and most of the heavy ones sink to the bottom; the 
end state is similar to the one seen in Figure \T2\ for the low- 
gravity simulations. Due to the strong shakes, the particles 
suffer large displacements, but, in a statistical sense, the two 
set of particles are segregated. A density segregation is then 
observed, although it is not complete. 

The results of the simulations under a low-gravity envi- 
ronment are presented in Figure [T5b . The lines for the light 
and heavy particles median height do cross for the three 
studied floor velocities {v_floor — {0.03,0.05,0.1} m/s), 
although for the lowest velocity the simulations do not last 
long enough to reach the stable stage where the median 
heights reach almost a stable value. 

Note that in both gravity environments, the density seg- 
regation is effective for floor's velocity over a threshold sim- 
ilar to the ones of the size segregation effect of Section [3] 



5 PARTICLE LIFTING AND EJECTION 

Let us consider the following simple experiment: we have a 
layer of material that is uniformly shocked from the bot- 
tom. The motivation of this experiment is to consider what 
would happen if a seismic wave, generated somewhere in a 
body and propagating through it, reaches another region of 
the body from below. What would happen with material 
deposit on the surface? Let us take into account three dif- 
ferent materials: a solid block, a compressible fluid and a set 
of grains. The outcome of the experiment will be different 
depending on the material. When the seismic wave knocks 
the solid block, the block is pushed upward. It starts to move 
upward, forming a gap between the layer's bottom and the 
floor. In the case of a layer of compressible fluid, an elas- 
tic p-wave is transmitted through it, producing compression 
and rarefaction of the material. 

But, what happens in the case of a layer of grains? Be- 
fore presenting the results of some simulations, let us re- 



consider the simulations of Newton's cradle with Hertzian 
viscoelastic spheres. We have seen that after the flrst par- 
ticle knocks the second one from the right, all the particles 
move to the left. Particle #4, the last one on the row, moves 
faster, the next one to the right moves slower and so forth. 
Therefore, the whole set of particles move in the same direc- 
tion, but they do not do it as a compact set, the particles 
separate from each other. 

We perform a first set of simulations with a homoge- 
neous set of particles. A 3D box with a base of 7.5 x 7.5 m 
is filled with 15 x 15 = 225 particles glued to the bottom, 
with a radius R — 0.25 m. We create 2744 particles with a 
mean radius Ri — 0.25 m, standard deviation a = 0.01 m 
and density p = 3000 kg/m? . To generate the initial condi- 
tions for the simulations, these particles are located a few 
cm from the bottom and they free fall under the different 
gravity environments until they settle down. 

Elastic parameters of all the particles are similar to the 
ones used in the previous tests for all the particles: Y = 
10^" Pa, A = 10"^ s"\ V = 0.3, k = 0.4, = 0.6, K = 
10^ Pa. 

With the initial conditions generated above, we run the 
following experiment: after a given time (taep), the fioor is 
vertically displaced at a certain speed {v floor) for a short 
interval {dt shake), only one time. Two gravity environments 
are used for the simulations: Earth's surface and the low- 
gravity environment of Itokawa. For the Earth's simulations 
we use the following set of parameters: tsep = 1 s, Vfioor ~ 
{1,3,10} m/s, dtshake = 0.1 s. At every snapshot, we sort 
the particles by their height respect to the fioor, and we 
compute the height of the particles at the 10% (/iio) and 90% 
(^9o) percentile. In Figure[l4]we plot the difference of these 
two quantities (/190 — /iio) for the different ffoor velocities. 
We observe that these differences increase with time up to a 
certain instant when the particles fall back. Therefore, the 
particles are not moving as a compact set, rather, the upper 
particles are moving faster and the particles separate from 
each other. The upper particles can reach velocities larger 
than the floor's velocity; e.g. in the case of v floor = 10 m/s, 
the 10% fastest particles reach velocities of ~ 17 m/s just 
after the end of the floor's displacement. We observe that 
the upper particles are lifted at considerable heights before 
they fall back. 

Similar results are obtained in low-gravity simula- 
tions, using the following set of parameters: taep ~ 10 s, 
Vfioor = {0.01,0.03,0.1} m/s, dtshake = 0.1 s. The up- 
per particles move faster and they can reach velocities up 
to ~ 0.02,0.05,0.2 m/s with respect to the floor veloci- 
ties. Note that the escape velocity in this environment is 
Vesc ~ 0.17 m/s, therefore the fastest ejection velocities of 
the lifted particles are higher than Vesc- We run another 
experiment: on top of the layer of particles with mean ra- 
dius -Ri ~ 0.25 m, we deposit a layer of 2700 smaller parti- 
cles, with mean radius 7?2 = 0.1 m and standard deviation 
(J = 0.01 m. The rest of the physical parameters are the same 
as for the bigger particles. The aim of this experiment is to 
check whether the small particles are ejected with higher 
velocities than the big ones. As in the previous simulations, 
we order the particles in increasing height. We compute the 
height of the 90% percentile of the big (/if,, 90) and small par- 
ticles (/is,9o)- Although the small particles on top of the big 
ones tend to separate, the differences in the velocities are 
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Figure 14. Lifting of particles under Earth gravity. At every 
snapshot, we sort the particles by their height respect to the floor, 
and we compute the height of the 10% (/iio) and 90% (/igo) per- 
centile. We plot the difference of these two quantities (hgo — /iio) 
for the different floor velocities. 
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Figure 15. The maximum height of the particles as a function 
of the simulated time for the case of the low-gravity environment 
and different floor velocities. 



relatively small. There is no significant ejection of the small 
particles. 

Another relevant result regarding the lifting and ejec- 
tion of particles from the surface due to an incoming shock 
from below, can be obtained from the Brazil nut effect sim- 
ulations presented in Section [S] In the animations produced 
with a sequence of snapshots for the simulations where the 
segregation process was effective, we observe many particles 
lifted at considerable heights. In Figure [15] we plot the max- 
imum height of the particles as a function of the simulated 
time for the case of the low-gravity environment and dif- 
ferent floor velocities. Note that the ejection velocities the 
fastest particles can acquire are comparable to the floor's 
displacement velocities, and even, a little bit higher. For a 
floor velocity of 0.1 m/s, the particles can reach an ejection 
velocity higher than the escape velocity at the surface. 

Taking into consideration the previous results, we con- 
clude that a layer shocked from below would produce the lift- 
ing of particles at the surface if the displacement of the bot- 
tom exceeds a certain velocity threshold. Particles can ac- 
quire vertical velocities comparable to the displacement ve- 
locity of the bottom. For very low-gravity environments, this 



velocity could be comparable to the escape velocity at the 
surface. The particles could enter in sub-orbital or orbital 
flights, creating a cloud of gravitational weakly bounded par- 
ticles around the object. 



6 GLOBAL SHAKING DUE TO IMPACTS 
AND EXPLOSIONS 

In the previous sections we have shown that several physi- 
cal processes can occur in a layer of granular media when it 
is shocked from below: size and density segregation, lifting 
and ejection of particles. A big quake in a distant point could 
produce such a shock. The quake could be produced by an- 
other small object impacting the body or by the release of 
some internal stress. Interplanetary impacts typically occur 
at velocities of several km/s. These are hypervelocity im- 
pacts, i.e. impacts with velocities that are above the sound 
speed in the target material, which give rise to physical de- 
formation of the target, heating and shock waves spreading 
out from the impact point. The DEM algorithms described 
above can not successfully reproduce these set of phenom- 
ena. Therefore, we have to implement a different approach 
if we are interested in understanding the effect of an impact 
induced shock wave passing through a granular media. 

Let's consider a fcm-size agglomerated body, formed by 
many m-size boulders. We raise the following question: what 
happens if a small projectile impacts in such an object at 
distances far from the impact point? Or alternatively, what 
happens if a large amount of kinetic energy is released in 
a small volume close to the surface of such an object? In 
order to answer these questions we run the following set of 
simulations. We fill a sphere of radius 250 and 1000 m with 
small spheres of a given size range, using the configurations, 
number of moving particles, and total mass listed in Table 
[2] For each sphere, we fill the volume with two different dis- 
tributions of small spheres: one with ~ 90, 000 particles and 
another one with a larger number of particles ~ 700, 000. We 
try to make the total mass of the moving particles similar 
for each of the studied radii. 

A time step of dt = 10~* s is used in all the simulations. 
The simulations are run in a cluster with Intel Xeon multi- 
core processors (Model E5410, at 2.33 GHz, with 12MB 
Cache). For cases B and D we use up to 8 cores. In these 
simulation of 10 s takes ~ 20 hr of CPU-time in 
each core. 

Since we can not successfully simulate the physics of a 
hypervelocity impact during the very short initial stages, we 
implemented another approach. At a given point on the sur- 
face we select a certain number of particles of the body that 
are close to this place. Each particle has at the beginning 
of the simulation a velocity along the radial vector toward 
the centre. We substitute the impact by a near-surface un- 
derground explosion, where several particles are released at 
a given speed. For each set of configurations listed in Table 
(2] we run simulations with initial particle velocities of 100 
m/s and 500 m/s. These velocities are well below the sound 
speed in the target material. 

Since these initial conditions would correspond to a 
stage after the impact where some energy has already been 
spent in the compression, fracturing and heating of the tar- 
get material, we cannot equal the sum of the kinetic en- 
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Table 2. Parameters for the simulations of underground explosions 



Case 




A 


B 


C 


D 






R iiri 1 1 1 Q 


R P ri 1 11 1; 


R a H 1 1 1 Q 


R 1 1 iQ 






250 m 


250 m 




lono m 

A-\J\J\J ill 


Size range of spheres (m) 




2.5 - 12.5 


1 - 10 


10 - 50 


5 - 25 


Number of particles 




88570 


783552 


89144 


688443 


Porosity 




0.31 


0.22 


0.31 


0.31 


Total Mass (lO^^feg) 




0.135 


0.152 


8.66 


8.63 


Escape velocity at surface (m/s) 




0.269 


0.285 


1.075 


1.073 


Number of initially moving particles 




10 


140 


10 


200 


Mass of moving particles (lO^fcg) 




21 


21 


599 


608 


Energy-equivalent projectile radius (m) for v = 


lOOm/s 


9.88 


0.88 


2.67 


2.68 


Energy- equivalent projectile radius (m) for v = 


500m/s 


2.57 


2.56 


7.82 


7.85 


Momentum-equivalent projectile radius (m) for 


V = lOOm/s 


3.26 


3.23 


9.84 


9.89 


Momentum- equivalent projectile radius (m) for 


V = 500m/s 


5.53 


5.52 


16.83 


16.92 


Ratio Kinetic Energy / Potential Energy for v 


= 100 m/s 


36 


29 


1 


1 


Ratio Kinetic Energy / Potential Energy for v 


= 500 m/s 


907 


710 


25 


26 


Specific energy Q* (J/kg) for v = 100 m/s 




0.79 


0.69 


0.35 


0.35 


Specific energy Q* (J/kg) for v = 500 m/s 




20 


17 


8.6 


8.8 



ergy of the moving particles with the kinetic energy of the 
impactor. However, we can provide a lower limit to the ki- 
netic energy of the impactor by assuming efficiency factor 
ske = 1, or a corresponding lower limit of the impactor 
size for a given impact velocity. In Table (2] we also present 
the radius of the equivalent projectile for the two set of ini- 
tial particle velocities, assuming an energy efficiency factor 
ekb = 1 and an impact velocity of 5 km/ a. For lower values 
of the efficiency factor, the projectile radius would scale with 
e^^^. As we have seen in the simulations of the Newton's 
cradle with viscoelastic interactions, there is a considerable 
loss of kinetic energy after a series of collisions, although 
the total linear momentum is conserved. As far as we know, 
there is very limited data on the transfer of momentum in 
hypervelocity impacts, and we do not know the efficiency 
factor of this transfer (clm). A similar estimate of the lower 
limit for the impactor size can be done by assuming a mo- 
mentum efficiency factor eLM ~ 1 and an impact velocity of 
5 km/s. In Table [2J we present the radius of the equivalent 
projectile for the two set of initial particle velocities. The 
projectile radius would scale with e^^^ . 

The location of the explosion is always at the sur- 
face and with angular coordinates (latitude = 45 deg , 
longitude — 45 deg). In Figure [TS] we present snapshots 
showing the propagation of the wave into the interior, by 
using slices passing through the centre of the sphere, the 
explosion point and the poles. Figure [TU] a and b correspond 
to the simulations with body radius of 250 m, the largest 
number of particles (A'' — 783552) and particles velocities 
of 100 m/s (case B-lOO). Snapshot a is at 0.4 s after the 
explosion and b is at 2 s. The particles are coloured using a 
colour bar that scales with the modulus of the velocity. On 
the other hand. Figure 116b and d correspond to the simu- 
lations with body radius of 1000 m, the largest number of 
particles (A'^ = 688443) and particles velocities of 500 m/s 
(case D-500). Snapshot c is at 3 s after the explosion and d 
is at 6 s. In the supplementary material we present movies 
of these simulations, (movie? corresponds to the case B-lOO 




Figure 16. Snapshots of the sphere explosions simulations. These 
arc slices passing through the centre of the sphere, the explosion 
point and the poles. Snapshots a) and b) correspond to the simu- 
lations with body radius of 250 m, the largest number of particles 
(N = 783552) and particles velocities of 100 m/s (case B-lOO). 
Snapshot a is at 0.4 s after the explosion and b is at 2 s. The parti- 
cles are coloured using a colour bar that scales with the modulus 
of the velocity. Snapshots c) and d) correspond to the simula- 
tions with body radius of 1000 m, the largest number of particles 
(N = 688443) and particles velocities of 500 m/s (case D-500). 
Snapshot c is at 3 s after the explosion and d is at 6 s. (see movies 
in the supplementary material) 

m/s, movie8 to case B-500 m/s, movieO to case D-lOO m/s, 
and movielO to case D-500 m/s. In the movies we observed 
the variation of the velocity of the particles in a slice pass- 
ing through the centre of the sphere, the explosion point and 
the poles. The particles are coloured using a colour bar that 
scales with the modulus of the velocity. 

We note that a shock front with a spherical shape prop- 
agates to the interior from the explosion point. On the sur- 
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face, there appears a layer of fast moving particles that ex- 
tends until it intersects with the spherical front, creating in- 
side the volume limited by the surface layer and the spherical 
front, a cavity of slow moving particles. The velocity of the 
propagation front has a weak dependence on the velocity of 
the initial particles. For example, in the simulations of the 
smaller body (case B-lOO), the propagation shock requires 
1.8 s to reach the antipodes of the explosion point, implying 
a velocity of 278 m/s. In the case B-500, the required time 
is 1.2 s, and the velocity 416 m/s. For the largest body, the 
figures are: case D-lOO: time 9.6 s, velocity 208 m/s; case 
D-500: time 5.8 s, velocity 435 m/s. Although there is an 
increase in the initial velocity of the moving particles of a 
factor of 5 among the cases, the velocity of the propagation 
shock has an in increase of ~ 2. The velocity of the propaga- 
tion shock is quite constant while the shock travels through 
the interior. 

We are interested in the effects of the explosion at large 
distances from the explosion point. The body is divided in 8 
quadrants. The explosion occurs on the surface at the cen- 
tre of the first quadrant (in Cartesian coordinates the first 
quadrant is:2;>0&j/>0&2:>0; and the explosion 
point is a,t: X = y — z = R/\/i, R - radius). We analyse the 
distribution of ejection velocities of the particles close to the 
surface (r > 0.8i?) on the other 7 quadrants. Histograms of 
these distributions are presented in Figure [17] a and b. In 
Figure 117b there are two overlapping histograms which cor- 
respond to the cases B-lOO and B-500, while in Figure [iTb . 
they correspond to the cases D-lOO and D-500. A vertical 
line marking the escape velocity for each body is included 
in the plots. Note that for the smallest object and for both 
initial velocities, there is a significant fraction of particles 
that acquire ejection velocities over the escape limit. Con- 
sidering the total fraction of particles with velocities over 
this threshold (not only the ones near the surface), we ob- 
tain values of 18% in the case B-lOO, and 81% in the case 
B-500. In the case of the largest body, there is a significant 
fraction of escaping particles only for the largest initial ve- 
locity. The total fraction of escaping particles are 0.6% in 
the case D-lOO, and 100% in the case D-500. For the sim- 
ulations with initial velocities of 500 m/s, there is a total 
disruption of both bodies (> 50% of the mass is ejected 
at velocities over the escape one). It is out of the scope of 
this paper to derive the disruption laws for this type of ex- 
periments; we just mention that with a set of experiments 
like the previous ones, we could obtain the kinetic energy 
threshold over which the explosions lead to a total disrup- 
tion of the body as a function of size. In Table [5] we also 
include the ratio between the kinetic energy of the moving 
particles over the potential energy of the body and the spe- 
cific energy (defined as the de posited energy per unit mass) . 
iHousen and Holsappld l|l990h have defined the critical spe- 
cific energy (Q*) as the energy per unit ma ss necessary to 
catastrophically disrupt a body. iRvanI (|2000t ) presents a plot 
comparing different estimates of Q* by several authors as a 
function of the target radius. Let us note the fact that the 
largest body {R — 1000 m) is more disrupted than the small- 
est body {R = 250 m) , although the specific energy is lower, 
it is in agreement with the dip in the Q* vs R plot (|Rvanl 
l|2000l )') in this radius range. 

Except for the case of low velocity explosions for the 
large body, in all the other simulations, a fraction of the 



a b 




Velocity (m/s) Velocity (m/s) 

Figure 17. The distribution of the ejection velocity of the parti- 
cles for the simulated explosion, a) Simulations with body radius 
of 250 m and the largest number of particles (A^ = 783552) (case 
B). Two histograms are presented for initial velocities of 100 and 
500 m/s. b) Simulations with body radius of 1000 m and the 
largest number of particles {N = 688443) (case D). Two his- 
tograms for initial velocities of 100 and 500 m/s. In each plot, a 
vertical dashed line is drawn at the value of the escape velocity 
at the surface. 

near surface particles far from the explosion point acquire 
velocities over the escape one (see Figure [TT]). Therefore, an 
explosion would induce the ejection of particles from the sur- 
face at low velocities. These particles could either enter into 
orbit around the body or slowly escape from it, producing a 
cloud of fine particles that may take many days before disap- 
pearing. This result is complementary to the one obtained 
in Section [5] regarding the lifting and ejection of particles 
produced by a shake coming from below the surface. 

In the case of the smallest body, even the low velocity 
explosions would induce displacement velocities over several 
tenths of m/s on many near surface particles far from the 
explosion point (see Figure I17p . This displacement would 
produce a shake coming from below, similar to the shakes 
simulated in Section O The surface gravity of the smallest 
body is similar to the surface gravity used in the low-gravity 
simulations of Section [S] and for the largest body the con- 
ditions are similar to the simulation of Eros. Looking back 
to Figure [161 we conclude that explosion events like the one 
produced in our simulations would be enough to induce the 
shaking required to produce size and density segregation on 
the surface of these bodies. 

This process of shaking the entire object after an im- 
pact is suitable for small bodies where the escape velocity is 
comparable to the impact induced displacement velocity at 
large distance from the impact point. Further work should 
study up to which body sizes the shaking process is expected 
to occur. 



7 CONCLUSIONS AND APPLICATIONS OF 
THE RESULTS 

The main objective of this paper is to present the applica- 
tions of Discrete Element Methods for the study of the phys- 
ical evolution of agglomerates of rocks under low-gravity en- 
vironments. We have presented some initial results regard- 
ing process like size and density segregation due to repeated 
shakings or knocks, the lifting and ejections of particles from 
the surface due to an incoming shock and the effect of a sur- 
face explosion on a spherical agglomerated body. We recall 
that our shaking process is due to repeated set of knocks. 
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The main conclusions of these preliminary results are: 

• A shaking induced size segregation -the so-called Brazil 
nut effect-does occur even in the low-gravity environments 
of the surface of small Solar Systems bodies, like fem-size 
asteroids and comets. 

• A shaking induced density segregation is also observed 
in these environments, although it is not complete. 

• A particle layer shocked from below would produce the 
lifting of particles at the surface, which can acquire vertical 
velocities comparable to the surface escape velocity in very 
low- gravity environments. 

• A surface explosion, like the one produced by an impact 
or the release of energy by the liberation of internal stresses 
or by the re accommodation of material, would induce a 
shock transmitted through the entire body, and the ejection 
of surface particles at low velocities at distances far from 
the explosion point. This process is only suitable for small 
bodies. 

The application of these results to real cases will be the 
subject of further papers, but we foresee some situations 
where the results presented here will be relevant: 

• The internal structure of asteroid Itokawa and similar 
small asteroids formed as an agglomerate of m-size parti- 
cles, and the relevance of the Brazil nut effect produced by 
repeated impacts. 

• The non-uniform distribution of active zones in comets, 
like P/Hartley 2, and the internal density segregation of icy 
and rocky boulders produced by shakes caused by explosions 
and impacts. 

• The formation of dust clouds at low escaping velocities 
after an impact onto a fcm-size asteroid. 

The supplement online material can be accesed at: 
http : //www. astr onomia . edu . uy /Publi cat ions/ ^ 
|Taiicredi /Granular_Physi cs/ 1 
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SUPPLEMENTARY ONLINE MATERIAL FOR 
"GRANULAR PHYSICS IN LOW-GRAVITY 
ENVIRONMENTS USING DEM" 

Hereby you will find a set of movies included in the article 
"Granular physics in low-gravity environments using DEM" 
by Tancredi et al. (MNRAS, 2011). 

The supplement online material can be accesed at: 
http : / /www ■ astronomia . edu . uy/Publ i cat i ons/ 1 
Tancredi/Granular_Physics/ 



SIZE SEGREGATION (THE BRAZIL NUT 
EFFECT) SIMULATIONS 

A 3D box is constructed with elastic mesh walls. The box 
has a base of 6 x 6 m and a height of 150 m. A set of 12 x 12 m 
small balls of radius _Ri = 0.25 m are glued to the floor. The 
big ball has a radius R2 = 0.75 m, and on top of it, there 
are 1000 small balls with radii _Ri ~ 0.25 m. 

The floor is displaced with a staircase function as de- 
scribed in the paper with different velocities. 

We present movies for two set of simulations: a) under 
Earth's gravity (surface gravity g = 9.81 m/ s^) and a floor's 
velocity {v floor ~ 5 m/s), b) in a low-gravity environment 
[g = 10~* vn/s^) and {vfioor ~ 0.05 m/s). 

moviel.avi is a movie with all the spheres drawn and 
movie2.avi with the small spheres erased for the first simu- 
lation. The movies correspond to 100 seconds of simulated 
time and 50 shakes. 

While movie3.avi and movie4.avi correspond to the sec- 
ond one. The movies correspond to 10000 seconds of simu- 
lated time and 667 shakes. 



DENSITY SEGREGATION SIMULATIONS 

A 3D box similar to the previous one is created, with a 
6 X 6 m base and a height of 150 m. The box is con- 
structed with elastic mesh walls. On the floor we glue a 
set of 12 X 12 small balls of radius Ri — 0.25 m and den- 
sity p = 2000 kg/m^. There are 500 light balls with radii 
Ri ~ 0.25 m and density p = 500 kg/m^. On top of them, 
there are 500 heavy balls with similar radii and density 
p = 2000 kg/m^. At the be ginning of the simulations the 
balls are placed sparsely, the light balls at the bottom and 
the heavy ones on top. They free fall and settle down before 
starting the floor shaking. 

The floor is displaced with a staircase function in a sim- 
ilar way as in the previous set of simulations. 

We present movies for two set of simulations: a) under 
Earth's gravity (surface gravity g — 9.81 m/ s^) and a floor's 
velocity {vfioor = 3 m/s), b) in a low-gravity environment 
[g = 10~* m/s^) and {vjioor = 0.05 m/s). 

movie5.avi is a movie of the flrst simulation, while 
movie6.avi corresponds to the second one. The movie5 cor- 
responds to 1000 seconds of simulated time and 500 shakes, 
while the movie6 corresponds to 20000 seconds of simulated 
time and 1333 shakes. 

Note that in these movies the camera moves with the 
floor, therefore it seems that the floor is always located in 
the same position, but it really is moving with the staircase 
function described above. 



In movie5 the density segregation is not reached; while 
in movie6, most of the light particles move to the top and 
most of the heavy ones sink to the bottom. 



GLOBAL SHAKING DUE TO IMPACTS AND 
EXPLOSIONS 

We consider a km-size agglomerated body, formed by many 
small size boulders. We fill a sphere of radius 250 and 1000 
m with ~ 700,000 small spheres of a given size range (1- 
10 m-size boulders in the case of the small body, and 5-25 
m-size boulders for the big body). 

At a given point on the surface we select a certain num- 
ber of particles of the body that are close to this place. 
Each particle has at the beginning of the simulation a veloc- 
ity along the radial vector toward the centre. The location 
of the explosion is always at the surface and with angular 
coordinates {latitude = 45deg , longitude — 45deg). We 
run simulations with initial particle velocities of 100 m/s 
and 500 m/s. 

In the movies we present snapshots showing the propa- 
gation of the wave into the interior. These are slices passing 
through the centre of the sphere, the explosion point and 
the poles. The particles are coloured using a colour bar that 
scales with the modulus of the velocity. 

movie7.avi and movie8.avi correspond to the simulation 
with a body of radius 250 m, N = 783552 small particles and 
140 particles with initial velocities of 100 m/s (case B-lOO) 
and 500 m/s (case B-500), respectively. 

movie9.avi and moviel0.avi correspond to the simula- 
tion with a body of radius 1000 m, N = 688443 small parti- 
cles and 200 particles with velocities of 100 m/s (case D-lOO) 
and 500 m/s (case D-500). All the movies correspond to 10 
seconds of simulated time. 



